function [xc,resid,normR,dumped] = dampedNewton(fhandle,x0,options,exDat0,exDat1)

tol = 1e-8;
relax = .2;
maxIts = 1e2;
nExDat = 0;
relTol = 1e-9; 
if nargin > 2
    if isfield(options,'tol')
        tol = options.tol; % stop tolerance on residual
    end
    
    if isfield(options,'relTol')
        relTol = options.relTol;
    end
    
    if isfield(options,'relax')
        relax = options.relax; % newton relaxation/damping factor
    end
    
    if isfield(options,'maxIts')
        maxIts = options.maxIts;
    end
    
    nExDat = nargin - 3;
end

opts.a = 0;

switch nExDat
    case 0
        [r,j] = feval(fhandle,x0);
    case 1
        [r,j] = feval(fhandle,x0,exDat0);
    case 2
        [r,j] = feval(fhandle,x0,exDat0,exDat1);
        j = fdJac(fhandle,x0,opts,exDat0,exDat1);
%         norm(j - j0)
end

xc = x0;

% xc = x0 - relax*(j0\r0);

done = 0;
its = 1;
normR = zeros(1,maxIts);
normR(1) = norm(r);

while ~done
    its = its + 1;
    xc = xc - relax*(j\r);
    
    switch nExDat
        case 0
            [r,j] = feval(fhandle,xc);
        case 1
            [r,j] = feval(fhandle,xc,exDat0);
        case 2
            [r,j] = feval(fhandle,xc,exDat0,exDat1);
            j = fdJac(fhandle,xc,opts,exDat0,exDat1);
%             norm(j-j0)
    end
    
%     norm(abs(r)) < tol
%     
%     its > maxIts
    
    done = (norm(abs(r)) < tol) + (its > maxIts) + ...
        (norm(r)/normR(1) < relTol);
    normR(its) = norm(r);
    
%     fprintf('\nNorm(R) = %f\n',norm(r))
end

% if (norm(r)/normR(1) < relTol)
%     fprintf('\n Rel Tol reached with %e\n',norm(r)/normR(1));
% end
% 
% if (norm(abs(r)) < tol)
%     fprintf('\n Abs Tol reached with %e\n',norm(r));
% end

% fprintf('\nNewton Finished in %d iterations\n',its);
dumped = 0;
if its > maxIts
    disp('Newton reached max iterations and dumped out')
    dumped = 1;
end

resid = r;
    

    
